#!/usr/bin/perl -w
# -----------------------------------------------------------------------+
#                                                                        |
# GET CLONE PAIRS                                                        |
#                                                                        |
# -----------------------------------------------------------------------+
# get_clone_pairs.pl                                                     |
#                                                                        |
# AUTHOR: James C. Estill                                                |
# CONTACT: jestill@uga.edu                                               |
# STARTED: 3/11/06                                                       |
# UPDATED: 3/13/06                                                       |
# DESCRIPTION:                                                           |
# Generate a single record per clone for a set of clones that have been  |
# seqeunced from both ends.The resulting "assembled" clones are reverse  |
# complementd on the reverse read, stored in a MySQL database, written   |
# to a new fasta file, and formatted for nucleotide BLAST.               |
# The input file can be any sequence data format that is recognized by   |
# BioPERL.                                                               |
#                                                                        |
# DEPENDENCIES:                                                          |
#  - BioPerl                                                             |
#  - DBI                                                                 |
#  - GetOpt                                                              |
#                                                                        |
# USAGE:                                                                 |
#  - get_clone_pairs.pl -i infile -f infile_format -u username           |
#                       -p password -t table_name [-q]                   |
#                                                                        |
# OTHER VARIABLES TO CONSIDER PASSING TO THE COMMAND LINE:               |
# These would not be required and would have default values              |
#  - l Default length of Ns      || "100"                                |
#  - n Database Name             || "maize_hmpr_msll"                    |
#  - x TextName Delim            || "."                                  |
#  - w Forward String Identifier || "f"                                  |
#  - v Reverse String Identifier || "w"                                  |
#                                                                        |
#------------------------------------------------------------------------+

# -------------------------+
# INCDLUDES                |
# -------------------------+
use Bio::SeqIO;            # Reading and writing of seq files in different formats
use Bio::Seq;              # Seq objects
use DBI();                 # Database interface
use Getopt::Std;           # Get options from the command line

# --------------------------+
# VARIABLES                 |
# --------------------------+
my $SeqNum = 1;
my $gap = "NNNNNNNNNN";
my $BlastDir = "/home/jestill/blast/msll_hmpr/"; # Directory for the formatted blast database

# --------------------------+
# COMMAND LINE VARIABLES    |
# --------------------------+
my %Options;                 # Hash to hold the options from the command line
getopts('i:f:u:p:t:q', 
	\%Options);
my $quiet;
$infile = $Options{i};        # Full infile path to work with
$infileformat = $Options{f};  # Format of the infile
                              # Any valid bioperl format will work
$DbUserName = $Options{u};    # Database user name
$DbUserPassword = $Options{p};# Database password
$TableName = $Options{t};     # Name of the table that will be created 
$quiet = $Options{q};         # Load the quiet variable if defined at the command line


# If the required variables are not given at the command line 
# then give the user a warning with usage and exit
if (! $infile || ! $infileformat || !$DbUserName  || ! $DbUserPassword )
{
    print "\nUsage:  - get_clone_pairs.pl -i infile -f infile_format -u username\n";
    print "                               -p password -t table_name [-q]\n";
    print "\nUsage: GenBankSeqUpdate -u user -p password -d data table list\n";
    print " -u REQUIERD username for database connection\n";
    print " -p REQUIRED password for database connection\n";
    print " -i REQUIRED infile\n";
    print "    The infile containing the sequence data to assemble\n";
    print " -f REQUIRED infile format\n";
    print "    The file format of the file to be assembled.\n";
    print "    (abi, ace, gcg, genbank, fasta, swiss, tigr etc.)\n";
    print " -t REQUIRED table name\n";
    print " -q OPTIONAL quiet flag\n";
    print "    In quiet mode no output is printed to the STDOUT\n";
    exit;
}


# -------------------------+
# GET FILE NAMES AND SET   |
# PATHS FOR DATA OUTPUT    |
# -------------------------+
my @file_name_parts = split( /\//, $infile );
my $LenArray = $#file_name_part;
my $process_file_name = $file_name_parts[$LenArray];
my $outfile = "$infile"."_pairs.fasta";

# -------------------------+
# SEQUENCE IO/DB HANDLES   |
# -------------------------+

# The sequence input file
# May need to specify here that thits is DNA
my $seq_in = Bio::SeqIO->new('-file' => "<$infile",
			     '-format' => $infileformat);

# CONNECT TO THE DATABASE
my $dbh = DBI->connect("DBI:mysql:database=maize_hmpr_msll;host=localhost",
		       $DbUserName, $DbUserPassword,
		       {'RaiseError' => 1});

# -------------------------+
# GENERATE TABLE FOR       |
# STORING OUTPUT DATA      |
# -------------------------+
# Test this as MSLL data. Will need to change to
# variable name based on the file name

$CreateTable = " CREATE TABLE ".$TableName."(".
    " rownum INT(10) NOT NULL AUTO_INCREMENT,".
    " clone CHAR(50),".
    " forseq LONGTEXT,".
    " revseq LONGTEXT,".
    " PRIMARY KEY (clone), ".
    " KEY (rownum))";
    
if (&does_table_exist($dbh, $TableName))
{
    # May want to warn that we are deleting a table here
    $dbh->do("DROP TABLE $TableName");
}

$dbh->do($CreateTable);

if (! $quiet) {print "The table ".$TableName." has been created\n";}


# -----------------------------------+
# APPEND CLONE DATA AND SEQUENCE     |
# DATA TO THE TABLE                  |
# -----------------------------------+
# This uses the reverse complement for the reverse sequence
# to put the forward and reverse sequences in the same frame
# It would be more effecient to run this as and INSERT/IGNORE
# statement in MySQL
while (my $inseq = $seq_in->next_seq())
{

    if (! $quiet) {print "Processing $SeqNum : ".$inseq->primary_id."\n";}
    my $info = $inseq->primary_id;
    @info_split = split(/\./, $info);
    my $clone_id = $info_split[0];
    my $seq_direction = $info_split[1];

    # ---------------------------+
    # INSERT CLONE ID TO RECORD  |
    # ---------------------------+ 
    # Append the clone name to the table tblMSLL_pairs
    # USING INSERT INGORE will ignore duplicate key problems
    my $AppendRecord = "INSERT IGNORE INTO ".$TableName.
	" SET clone='".$clone_id."'";
    $dbh->do($AppendRecord);

    # --------------------------+
    # GET REVERSE COMPLEMENT    |
    # FOR r READS AND UPDATE    |
    # SEQ DATA TO TABLE         |
    # --------------------------+
    if ( $seq_direction =~ "r" )
 
    {
	# Edited the msll file to get rid of ZMMBLd0018I04.r
	$UploadRevSeq = "UPDATE ".$TableName.
	    " SET revseq='".$inseq->revcom->seq."'".
	    " WHERE clone ='".$clone_id."'";
	$dbh->do($UploadRevSeq);
    }
    elsif ( $seq_direction =~ "f" )
    {
	$UploadForSeq = "UPDATE ".$TableName.
	    " SET forseq='".$inseq->seq."'".
	    " WHERE clone ='".$clone_id."'";
	$dbh->do($UploadForSeq);
    }
    else
    {
	print "An unexpected direction has been seen \n";
    }
    
#    # Place to kick out if this is a test run
#    if ($SeqNum == 100)
#    {
#	exit;
#    }
#
    $SeqNum++;

}

# -----------------------------------+
# EXPORT THE DATA FROM THE TABLE     |
# TO A NEW PAIRED FASTA FILE         |
# -----------------------------------+
# May want to create a separate FASTA file
# for unpaired clones

if (! $quiet) {print "Printing out to file $outfile\n";}
open (FILEOUT, ">".$outfile);

$NumRecords = &how_many_records($dbh,$TableName) ;

for ($loop_index = 1; $loop_index <= $NumRecords; $loop_index++) 
{

#    # PLACE TO EXIT FOR A TEST RUN
#    if ($loop_index == 100) {exit;}

    $Query = "SELECT * FROM $TableName WHERE rownum=$loop_index" ;
    $sth = $dbh->prepare($Query);
    $sth->execute() || print "\nERROR: Can not execute select for $TableName with rownum: $loop_index \n";
    
    while(@row_ref = $sth->fetchrow_array)
    {
	print FILEOUT ">".$row_ref[1]."\n".$row_ref[2]."$gap".$row_ref[3]."\n";
	if (! $quiet) {print "$row_ref[1]\n";}
    } # END OF THE WHILE LOOP TO GET THE DATASET FROM THE QUERY
    
} # END OF LOOP THAT RUNS DOWN THE ROWS as rownum with $loop_index


# -----------------------------------+
# GENERATE A BLAST DATBASE FOR THE   |
# PAIRED FASTA FILE                  |
# -----------------------------------+

my $CopyCommand = "cp $outfile $BlastDir$process_file_name"."_paired";
my $BlastCommand = "formatdb -i $BlastDir$process_file_name"."_paired".
    " -p F -t $process_file_name"."_paired -n $process_file_name"."_paired";

if (! $quiet) {print "Copying file to BLAST directory.\n";}
system ( $CopyCommand );

if (! $quiet) {print "Formating the BLAST database.\n";}
system ( $BlastCommand );

# -----------------------------------+
# ALERT THE USER THAT ALL IS DONE    |
# -----------------------------------+

print "PROCESSED FILE: $process_file_name\n";
print "\t  TABLE: $TableName\n";
print "\t SOURCE: $infile\n";
print "\t    OUT: $outfile\n";
print "\tBLASTDB: $BlastDir$process_file_name"."_paired\n";
print "$CopyCommand\n";
print "$BlastCommand\n";

exit;

# --------------------------
# SUBROUTINE TO CHECK TO SEE
# IF THE TABLE ALREADY EXITS
# --------------------------

sub does_table_exist
{
    my ($dbh,$whichtable) = @_;
    my ($table,@alltables,$found);
    @alltables = $dbh->tables();
    $found = 0;
    foreach $table (@alltables) 
    {
    	$found=1 if ($table eq "`".$whichtable."`");
    }
    # return true if table was found
    # return false if not found
    return $found;
}


# ------------------------------------+
# SUBROUTINE TO SEE HOW MANY RECORDS  |
# EXIST IN THE MYSQL TABLE            |
# ------------------------------------+
# CODE FROM
# http://lena.franken.de/perl_hier/databases.html

sub how_many_records
{

  my ($dbh,$whichtable) = @_;
  my ($result,$cur,$sql,@row);
  
  # trace what is going on for error checking purposes
  # DBI->trace( 3 );
  $sql = "select count(*) from $whichtable"; # run the query that counts the number of rows
  $cur = $dbh->prepare($sql);                # prepare the query as cur
  $cur->execute();                           # execute the query
  @row=$cur->fetchrow;                       # fetch the rows in the row array
  $result=$row[0];                           # the string result is the first row in the array
  # Code below added by JCE to see the number of rows that the program sees in the loop for debug purposes
  # print "\n=\n===\n===== Result in loop is $result \n=====\n===\n=\n";
  $cur->finish();                            # this finises the cur string ... kills it
  return $result;                            # return the result to the user
  
}
